Aperiodic-Order- Assisted Optical Nonlocality in Multilayered Hyperbolic 

Metamaterials 
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We show that hyperbolic electromagnetic metamaterials implemented as multilayers based on 
two material constituents arranged according to Thue-Morse aperiodic sequence may exhibit strong 
nonlocal effects, manifested as the appearance of additional extraordinary waves which are not 
predicted by standard effective-medium-theory (local) models. These effects can be associated with 
stationary points of the transfer-matrix trace, and can be effectively parameterized via the trace-map 
formalism. Interestingly, for certain parameter configurations, at a given wavelength and for two 
given material layers, these effects may disappear when the same layers are arranged in a standard 
periodic fashion, thereby providing a striking evidence of the major role played by aperiodic order. 
Our findings indicate that the (aperiodic) positional order of the layers constitutes an effective and 
technologically inexpensive additional degree of freedom in the engineering of optical nonlocality. 

PACS numbers: 42.25.Bs, 78.67.Pt, 78.20.Ci, 61.44.Br 



Electromagnetic (EM) metamaterials are artificial ma- 
terials composed of subwavelength dielectric and/or 
metallic inclusions in a host medium, which have at- 
tracted considerable scientific and applicative attention 
due to possibility to engineer anomalous properties (e.g., 
negative refraction) that are not observable in natu- 
ral materials [lj. Of particular interest are the so- 
called "hyperbolic" metamaterials 0, HI, characterized 
by nonmagnetic, uniaxally-anisotropic constitutive rela- 
tionships with both positive and negative components 
of the permittivity tensor. This yields a hyperbolic (as 
opposed to spherical, in conventional isotropic media) 
dispersion relationship, which allows for propagation of 
(otherwise evanescent) waves with large wavevectors, re- 
sulting in a high (theoretically unbounded) photonic den- 
sity of states. The reader is referred to ji- 12| for a sparse 
sampling of applications, ranging from nanoimaging to 
quantum nanophotonics and thermal emission. 

In what follows, we focus on multilayered hyperbolic 
metamaterials Q , implemented via stacking of alternat- 
ing subwavelength layers with negative and positive per- 
mittivities (e.g., metallic and dielectric, at optical wave- 
lengths). For this class, the effective medium theory 
(EMT) provides a particularly simple model in terms of a 
homogeneous, uniaxially-anisotropic permittivity tensor 
with components given by the Maxwell-Garnett mi xing 
formulas [13| . However, a series of recent papers [Fa - tlg 
have pointed out the limitations of this model in predict- 
ing nonlocal effects that can take place (even in the pres- 
ence of deep subwavelength layers) due to the coupling 
of surface plasmon polaritons (SPPs) propagating along 
the interfaces separating layers with opposite-signed per- 
mittivities. This may result, for instance, in the mispre- 
diction of additional extraordinary waves HI, 3] as well 
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as of the broadband Purcell effect 

We point out that typical multilayered hyperbolic 
metamaterials are based on periodic arrangements of the 
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FIG. 1. (Color online) Problem schematic, illustrating the 
2-D propagation of TM-polarized EM fields (with y-directed 
magnetic field) in the ThM-based hyperbolic metamaterial of 
interest (details in the text). 



layers [7]. In fact, the EMT model describing the lo- 
cal response is independent of the positional order of the 
layers, and depends only on the permittivities of the two 
constituents and their filling fractions 13[. However, one 
would intuitively expect the positional order of the layers 
to sensibly affect the nonlocal response. It seems there- 
fore suggestive to investigate possible nonlocal effects in- 
duced by aperiodic order inspired by the "quasicrystal" 
concept in solid-state physics [19| . 

Accordingly, in this Letter, we study the nonlo- 
cal response of hyperbolic metamaterials based on 
aperiodically-ordered multilayer superlattices. With spe- 
cific reference to the Thue-Morse (ThM) geometry [20j, 
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we identify certain nonlocal effects (in terms of additional 
extraordinary waves) that can only be attributed to the 
specific positional order of the material layers. 

Referring to the two-dimensional (2-D) y— independent 
geometry in Fig. [TJ we consider a multilayer superlattice 
obtained via the infinite repetition along the z— axis of 
a supercell composed of layers of two nonmagnetic, ma- 
terial constituents labeled with the letters "a" and "6" 
(with relative permittivities e a and ej,, and thicknesses 
d a and d , respectively), arranged according to the ThM 
sequence. Assuming as an initiator the sequence "qfr," 
this amounts to iterating the following inflation rules [20( 



a — > ab, b — > ba. 



(1) 



as shown schematically in the inset of Fig. Q] for the 
first iteration-orders n. In what follows, we study the 
time-harmonic [exp(— iut)] propagation of transversely- 
magnetic (TM) polarized EM fields as the iteration-order 
n increases so as to approach the aperiodic-order regime. 
For simplicity, we neglect material losses, as previous 
studies [18| have shown that they only mildly affect op- 
tical nonlocality. 

It is readily recognized that the first two iterations 
(n = 1,2) correspond to standard periodic multilayers 
(with period d = d a +db and 2d, respectively). Although, 
strictly speaking, the geometry in Fig. [T] is inherently 
periodic for any finite iteration-order n, throughout, we 
refer to the first two orders as the "standard periodic" 
cases. Moreover, we highlight that, given the structure 
of the inflation rule in (fTJ), any iteration-order of our ThM 
multilayer contains the same proportions of constituents 
a and b as the standard periodic case, and differs solely 
in the positional order of the layers. Accordingly, the 
Maxwell-Garnett mixing formulas for the parallel (|j, i.e., 
x,y) and orthogonal (_L, i.e., z) permittivity components 



x^y) 

m 



£ad a + £bdb 



e a 1 d a + e b 1 d b 



(2) 



yield the same EMT model for any iteration-order, which 
results in the dispersion relationship 
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where k x and k z indicate the x— and z— components, 
respectively, of the wavevector k (cf. Fig. [1]), and 
ko = uj/cq = 2-k/\q indicates the vacuum wavenumber 
(with Co and Ao denoting the corresponding wavespeed 
and wavelength). By suitably choosing the parameters in 
the mixing rules ([2]) so that £||£_i_ < 0, the dispersion rela- 
tionship in ([3]) , interpreted in terms of equi-frequency con- 
tours (EFCs) , assumes the anticipated hyperbolic charac- 
ter. Since the local EMT model in @ and © is identi- 
cal for any iteration-order, any possible difference in the 



(nonlocal) responses should solely be attributed to the 
different positional order of the material layers. 

Multilayers based on the ThM geometry have been 
widely studied in the past, in the form of dielec- 
tric/semiconductor photonic quasicrystals (see, e.g., [2ll — 
l2r3 | for a sparse sampling), and with main focus on 
the resonant-transmission, localization, omnidirectional- 
reflection, and bandgap properties. To the best of our 
knowledge, no previous attempt was made to study ThM- 
based hyperbolic metamaterials. Following a rather stan- 
dard approach (see [27} for more details), the exact dis- 
persion relationship pertaining to a nth-order ThM su- 
percell terminated by Bloch-type phase-shift walls (cf. 
Fig. [1]) can be compactly written as 



cos (k z D n ) 



Xn 



(4) 



where D n = 2 n ~ 1 d represents the total supercell thick- 
ness at the iteration-order n, and \ n denotes the trace 
(i.e., the sum of the diagonal elements) of the transfer- 
matrix that relates the tangential components of the EM 
fields at the supercell interfaces x = and x — D n (see 
[271 ] for more details). For the first two iterations n = 1,2, 
the trace can be straightforwardly calculated as 

X„ = 2cos(n(5 a )cos(n(5 b )- ( — + — ] sin(n<5 a )sin(n£&) , 

V76 7a/ 

(5) 

thereby recovering the familiar Bloch-type dispersion re- 
lationship of standard periodic multilayers (as in [T3 - 
[H), where 



Sa,b = kza,bd a ,b, la, b — 



£ a ,bko 



(6) 



with k za ,b = \/fco £ a,6 - k%, Im(k za ,b) > 0. For higher- 
order iterations, a particularly simple recursive calcula- 
tion procedure can be adopted, based on the trace-map 
[H [i| (see also ^ for details) 



Xn+2 = Xn (Xn+1 - 2) + 2, U > 1. 



(7) 



The exact dispersion relationship in ((4]) reduces to the 
local EMT model in Q in the limit d — >• 0, but it may 
significantly depart from that for finite (and yet subwave- 
length) layer thicknesses. In particular, we are interested 
in exploring possible nonlocal effects manifested as the 
appearance of additional extraordinary waves that are 
not predicted by the local EMT model in (|3]). From a 
mathematical viewpoint, this phenomenon is related to 
multiple (apart from sign) k x solutions of (j4j) for a given 
value of k z and u), which may occur if the trace \n is 
a nonmonotonic function of fcf. We are therefore led to 
study the stationary points of the trace Xn with respect 
to the argument k 2 . In particular, in order to better em- 
phasize the role of aperiodic order in the onset of these 
nonlocal phenomena, we deliberately focus on parame- 
ter configurations for which the hyperbolic metamaterial 
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FIG. 2. (Color online) (a) Transfer-matrix trace (blue-solid 
curve) as a function of k x and (b) corresponding EFC (within 
the first Brilluoin zone), for a hyperbolic metamaterial with 
e a = 6.83, Eb = —1.83, and d a = db = d/2 — 0.05Ao (i.e., 
e II = 2.5, £± — —5) at the n = 1 iteration-order (i.e., stan- 
dard periodic multilayer of period d) . Also shown (red-dashed 
curves) are the predictions from the local EMT model. Due 
to symmetry, only positive values of k x and k z are shown, (c), 
(d) Same as above, but for the n = 3 iteration-order (four Bril- 
louin zones shown for direct comparison). The dash-dotted 
lines highlight the correspondence between the zero of xi and 
the maximum of X3- 



arising from the first iteration n = 1 (i.e., a standard 
periodic multilayer) is well-described by the local EMT 
model in ([3]). This translates in the trace xi m © being 
well approximated by its second-order Taylor expansion 
in d, 
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Figure HJa) compares, for one such parameter con- 
figuration (given in the caption, and corresponding to 
£|| = 2.5 and e± = —5), the exact trace xi [cf- ©] and 
its quadratic approximation ([8]), showing a reasonable 
agreement. The corresponding exact [cf. (jH)] and local- 
EMT [cf. ©] EFCs are compared in Fig. C^b) within the 
first Brillouin zone < fc 2 < tt/c?. Again, we observe a 
good agreement (especially for smaller values of k z ) and, 
most important, a single branch, which yields a single 
mode that is propagating for k z > y/£\\ka ~ 0.32ir/d, 
and evanescent otherwise. We now move on to looking 
at higher-order iterations of the ThM geometry. From 
the trace-map in ([7]), we straightforwardly obtain 



with the overdot denoting differentiation with respect to 
the argument k 2 ,. This implies that the vanishing of the 
trace at a given iteration-order, i.e., Xn = is a sufficient 
condition for a stationary point Xn+2 — 0, and hence the 
presence of additional extraordinary waves, at a higher- 
order iteration. Therefore, if the trace xi admits a zero, 
then X3 should exhibit at least one stationary point. For 
the assumed parameter configuration, for which the local 
EMT model, and hence the quadratic approximation in 
©, holds reasonably well at the first iteration-order, the 
position k x o of such zero (and corresponding stationary 
point) admits a simple analytical estimate as 
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Xn+2 = Xn [XnX.n+1 + ^Xn (X 
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2)], n>l, (9) 
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which yields a real solution provided that s^k^d 2 < 2. 
In our case, this last condition is verified and, as can be 
observed from Fig. [UJa), the estimate in (JTUJ) is moder- 
ately accurate, yielding a 9% error with respect to the 
actual zero position k x0 — 0.4138d/7r calculated numer- 
ically. For the same parameter configuration, Fig. [2][c) 
shows the trace X3 a t the n — 3 iteration-order, from 
which a maximum at k x o can be observed. As a conse- 
quence, besides a branch that is still in good agreement 
with the local EMT prediction, the corresponding EFCs 
shown in Fig. [UJd) [within a spectral region covering four 
Brilluoin zones, in order to facilitate direct comparison 
with Fig. [Hh)] exhibit additional branches, resulting in 
two additional modes (extraordinary waves) that prop- 
agate for arbitrarily small values of k z , and degenerate 
into a single mode at the band-edges ki m ^ = rrm/(2d), 
m = 0,1,2,... 

Interestingly, we observe that the maximum corre- 
sponds to X3 = 2, which represents the band-edge con- 
dition. When substituted in © (together with X3 = 0), 
this implies that also %4 = at k x o- Moreover, we note 
from the trace-map ([7]) that this will also imply that 
Xn = 2, i.e., Xn — 0, for any n > 5. In other words, 
the additional extraordinary waves associated with the 
stationary point at k X Q will be retained by any arbitrar- 
ily high iteration-order of the ThM multilayer, and hence 
also in the limit for which the artificial periodicity en- 
forced by the Bloch-type phase-shift walls is washed out 
by the actual aperiodic order. 

Qualitatively similar considerations also hold for pa- 
rameter configurations characterized by ey < and 
e_L > 0, not discussed here for brevity (see [13] for de- 
tails). 

Incidentally, in previous works, the condition Xn — 
2 has also been associated with perfect transmission 
through finite-sized (along z) dielectric ThM multilayers 
sandwiched between two infinite (along z) homogeneous, 
isotropic media [21M23I [26l] . Similar to that scenario, also 
in our case the number of additional extraordinary waves 
grows exponentially with the iteration-order, and can be 
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FIG. 3. (Color online) (a) Numerically-computed field- 
magnitude (H y ) map in false-color scale, for a slab of hy- 
perbolic metamaterial (parameters as in Fig. [2} at the n = 1 
iteration-order (standard periodic multilayer of period d; 16 
unit cells shown) of thickness 0.75Ao, embedded in vacuum, 
and excited by a unit-amplitude, normally-incident (k z — 0) 
plane- wave. Thin grey lines delimit the slab and layer inter- 
faces, (b) Longitudinal cut (black-solid curve) at z = 0.05Ao, 
and exponential fit (red-dashed curve) of the evanescent de- 
cay inside the slab as predicted by the local EMT model [cf. 
JTT}]. (c) Transverse cut at x = 0.05Ao, with the material 
layers visualized with different colors/shades. 



determined by generalizing well-established approaches 
(see, e.g., [26|). This is, however, beyond the scope of 
our investigation here. 

Instead, we focus on an independent validation of our 
findings above. To this aim, for computational affordabil- 
ity, we study the TM plane-wave propagation through 
a slab of our ThM-based hyperbolic metamaterial (in- 
finitely long in the z— direction and 0.75Ao thick along 
x, with parameters as in Fig. [2]) immersed in vacuum, 
at various iteration-orders. Figure EJa) shows a field- 
magnitude map, numerically computed via a rigorous- 
coupled- wave-analysis (RCWA; see [27| for details), per- 
taining to the first iteration-order n = 1 (i.e., standard 
periodic multilayer) for normal incidence (k z =0). As 
can be observed also from the longitudinal (x) cut in Fig. 
E^b), the field is totally reflected, with only an evanes- 
cent decay inside the slab, which is accurately fitted by 
an exponential tail (red-dashed curve) with attenuation 
coefficient 
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FIG. 4. (Color online) As in Fig. but for the n = 3 
iteration-order (four unit-cells shown). Longitudinal and 
transverse magnetic cuts in (b) and (c) are at z — 0.05Ao 
and x = 0.375Aq, respectively. 



as predicted by the local EMT model in ([3]). The z-cut 
in Fig. [3tc) shows that the transverse field distribution is 
rather uniform (~ 10% variations) and weakly peaked at 
the centers of the layers. Overall, as expected, the local 
EMT model provides a satisfactory prediction. 

The response dramatically changes for the n = 3 it- 
eration order, as illustrated in Fig. [4] In this case, a 
standing-wave pattern is clearly visible inside the slab. 
Looking at the longitudinal cut in Fig. HJb), from the 
distance of two consecutive peaks (~ 0.241Ao) we can es- 
timate a propagation constant k x « OAlbd/ir which is in 
excellent agreement with the k x o = 0.4138d/7r value per- 
taining to the degenerate additional extraordinary wave 
predicted by the EFCs in Fig. HJd) for k z = (normal 
incidence). The a;— cut in Fig. Hie), markedly different 
from the standard-periodic-multilayer counterpart in Fig. 
[3Jc), shows a transverse field profile with much larger am- 
plitude variations, and with peaks at the interfaces be- 
tween positive- and negative-permittivity layers, thereby 
evidencing the nonlocal nature of this mode, due to the 
coupling of SPPs propagating along thee interfaces. 

Concerning other iteration-orders, not shown here for 
brevity (see (27J for details), we found that the n — 2 case 
(standard periodic multilayer of period 2c?) qualitatively 
resembles the response in Fig. O and is still well de- 
scribed by the local EMT model. Higher-order iterations 
analyzed (up to n — 5) are instead much more similar to 
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Fig. [4j clearly evidencing nonlocal phenomena, though 
with more complex standing-wave patterns due to the 
presence of several additional extraordinary waves. 

To sum up, we have shown that hyperbolic metama- 
terials implemented as multilayered based on the ThM 
sequence may exhibit strong optical nonlocality (man- 
ifested as the appearance of additional extraordinary 
waves) at any iteration-orders but the first two {n = 1, 2, 
corresponding to standard periodic multilayers). From 
a mathematical viewpoint, we associated these effects to 
stationary points of the transfer-matrix trace, and de- 
rived simple analytical design rules. We emphasize that 
different iteration-orders differ solely in the positional or- 
der of the constituent material layers. Moreover, the par- 
ticular structure of the ThM inflation rule in (JI]) ensures 
that, at any iteration-order, no more than two consecu- 
tive identical symbols may occur (e.g., aaa and bbb, or 
longer, sequences are forbidden) [20]. This implies that 
the nonlocal effects observed are not trivially attributable 
to an effective increase of the average layer thickness, but 
are genuinely assisted by aperiodic order. To the best 
of our knowledge, against the many implications and 
applications of aperiodic-order to optics and photonics 
(see, e.g., 3(1 3l| for recent reviews), this represents the 
first evidence in connection with optical nonlocality. Be- 
sides the inherent academic interest, from the application 
viewpoint, this constitutes an important, and technolog- 
ically inexpensive, additional degree of freedom in the 
engineering of optical nonlocality, which may be also be 
fruitfully exploited within the recently-introduced frame- 
work of nonlocal transformation optics [32]. We high- 
light that the ThM sequence was considered here only 
in view of its particularly simple inflation rule and as- 
sociated trace-map, which facilitate analytical treatment 
as well as direct comparison with standard periodic mul- 
tilayers, but the results are more general. In fact, one 
of the most intriguing follow-up study may be the sys- 
tematic design of deterministic aperiodic sequences, via 
suitable inflation rules and associated polynomial trace- 
maps [28|], yielding prescribed nonlocal effects. 
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